SGI Developer Toolbox 6.1/SGI Developer Toolbox 6.1 - Disc 4.iso/lib/mathlib/libconv/TRY/ornl_diir2d.f

Text File | 1994-08-02 | 5.2 KB | 168 lines

      subroutine ornl_diir2d(    f,incf,ldf,ifx0,n_fx,ify0,n_fy,
     $                g,incg,ldg,igx0,n_gx,igy0,n_gy,
     $                h,inch,ldh,ihx0,n_hx,ihy0,n_hy)
      double precision f(ldf,*)
      double precision g(ldg,*)
      double precision h(ldh,*)
      integer incf, ldf, ifx0, n_fx, ify0, n_fy
      integer incg, ldg, igx0, n_gx, igy0, n_gy
      integer inch, ldh, ihx0, n_hx, ihy0, n_hy
c-----------------------------------------------------------------------------
c
c   diir2d  performs a 2D convolution in the time domain :
c    h(i,j) = beta * h(i,j) + alpha * Sum.Sum[ f(k,l) * g(i-k,j-l) ]
c
c-----------------------------------------------------------------------------
c
c   PARAMETERS:
c
c    f:    Vector containing sequence "f"
c    incf:    Increment between two successive lines   of "f"
c    ldf:    Increment between two successive columns of "f"
c    ifx0:    Index of the first element of each column of "f"
c    n_fx:    Number of elements (lines) of each column of "f"
c    ify0:    Index of the first column of "f"
c    n_fy:    Index of the  last column of "f"
c
c    g:    Vector containing sequence "g"
c    incg:    Increment between two successive lines   of "g"
c    ldg:    Increment between two successive columns of "g"
c    igx0:    Index of the first element of each column of "g"
c    n_gx:    Number of elements (lines) of each column of "g"
c    igy0:    Index of the first column of "g"
c    n_gy:    Index of the  last column of "g"
c
c    h:    Vector containing sequence "h"
c    inch:    Increment between two successive lines   of "h"
c    ldh:    Increment between two successive columns of "h"
c    ihx0:    Index of the first element of each column of "h"
c    n_hx:    Number of elements (lines) of each column of "h"
c    ihy0:    Index of the first column of "h"
c    n_hy:    Index of the  last column of "h"
c
c-----------------------------------------------------------------------------
c
c   PLEASE NOTE:
c    Please note that the array pointers must all point to the first
c    element of the array "(ifx0,ify0)", "(igx0,igy0)" and "(ihx0,ihy0)"
c     If "f" for example is defined as
c        dimension f(-25:45,10:21)
c    Then "diir2d" must be called with the following parameters
c        call diir2d( f(-25,10),(45-(-25)+1),-25,45,10,21 ... )
c
c-----------------------------------------------------------------------------
c
c   This Fortran Subroutine written by
c    Jean-Pierre Panziera
c    Silicon Graphics Inc
c    September 27, 1991
c
c-----------------------------------------------------------------------------

      double precision zero, one
      parameter (zero = 0.0, one = 1.0)
      integer iy0, iy1, ify1, igy1, ihy1
      integer k1, k2, k, j, jf, jg, jh
      integer ix0, ix1, ifx1, igx1, ihx1, if0, ifx, ihx

c-----------------------------------------------------------------------------
    if( (n_hx .le. 0) .or. (n_hy .le. 0) ) return
c-----------------------------------------------------------------------------
c
c   Compute Y Boundaries
c
c-----------------------------------------------------------------------------
    ify1 = ify0 + n_fy - 1
    igy1 = igy0 + n_gy - 1
    ihy1 = ihy0 + n_hy - 1

    iy0 = min(max(ihy0, ify0-igy0), ihy1+1)
    iy1 = max(min(ihy1, ify1-igy0), ihy0-1)

    ifx1 = ifx0 + n_fx - 1
    igx1 = igx0 + n_gx - 1
    ihx1 = ihx0 + n_hx - 1

    ix0 = min( max(ihx0, ifx0-igx0), ihx1+1)
    ix1 = min(ihx1, ifx1-igx0)

c-----------------------------------------------------------------------------
c
c   Zero the Front columns
c
c-----------------------------------------------------------------------------
    do j = ihy0, iy0-1
        ihx = 1
        jh = 1 + (j-ihy0)
        do k = 1, n_hx
        h(ihx,jh) = zero
        ihx = ihx + inch
        end do
    end do
c-----------------------------------------------------------------------------
c
c   Init the Plain columns
c
c-----------------------------------------------------------------------------
    if0 = 1 + ((ix0+igx0)-ifx0) * incf
    do j = iy0, iy1
        ihx = 1
        jh = 1 + (j-ihy0)
        do k = ihx0, ix0-1
        h(ihx,jh) = zero
        ihx = ihx + inch
        end do
        ifx = if0
        jf = 1 + (j+igy0-ify0)
        do k = k, ix1
        h(ihx,jh) = f(ifx,jf)
        ihx = ihx + inch
        ifx = ifx + incf
        end do
        do k = k, ihx1
        h(ihx,jh) = zero
        ihx = ihx + inch
        end do
    end do
c-----------------------------------------------------------------------------
c
c   Zero the TAIL columns
c
c-----------------------------------------------------------------------------
    do j = iy1+1, ihy1
        ihx = 1
        jh = 1 + (j-ihy0)
        do k = 1, n_hx
        h(ihx,jh) = zero
        ihx = ihx + inch
        end do
    end do
c-----------------------------------------------------------------------------
    if( (n_gx .le. 0) .or. (n_gy .le. 0) ) return
c-----------------------------------------------------------------------------
c
c   Call the 1D Convolution
c
c-----------------------------------------------------------------------------
    do j = iy0, ihy1
        jh = j-ihy0 + 1
        k1 = max( iy0, j-(n_gy-1))
        k2 = j-1
        do k = k1, k2
        jf = k-ihy0 + 1
        jg = j-k + 1
        call dfir1d(    h(1,jf), inch, ihx0, n_hx,
     $                g(1,jg), incg,    0, n_gx,
     $                h(1,jh), inch, ihx0, n_hx,
     $                -one, one)
        end do
        call diir1d(    h(1,jh), inch, ihx0, n_hx,
     $                g(1,1 ), incg,    0, n_gx,
     $                h(1,jh), inch, ihx0, n_hx)
    end do
c-----------------------------------------------------------------------------
      return
      end